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O ' Abstract 

(N 

— ■ We present a variational solution of the Langevin field equation describing 

the nonequilibrium dynamics of a harmonically trapped Bose-Einstein con- 
densate. If the thermal cloud remains in equilibrium at all times, we find 
that the equation of motions for the parameters in our variational ansatz 
are equivalent to the Langevin equations describing the motion of a massive 
Brownian particle in an external potential. Moreover, these equations are 
coupled to a stochastic rate equation for the number of atoms in the con- 
densate. As applications of our approach, we have calculated the collisional 

t^ \ damping rates and frequencies of the low-lying collective excitations of a con- 

densate with repulsive interactions, and have obtained a description of the 
growth and subsequent collapse of a condensate with attractive interactions. 
We have found a good agreement with the available experimental results in 
both cases. 
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I. INTRODUCTION 

The experimental realization of Bose-Einstein condensation in dilute atomic gases |T] fj| . 
has led to a large increase in the amount of both experimental and theoretical research on 
these quantum systems. Various theoretical predictions regarding equilibrium and nonequi- 
librium properties of degenerate Bose gases, can now be directly compared with experi- 
mental data. Regarding the zero-temperature behaviour of a Bose-Einstein condensate, a 
great deal of the physics is well described by the Gross-Pitaevskii equation, i.e., a mean- 
field equation for the macroscopic wave function of the condensate. It has led to very good 
agreement with experimental results on, for example, the condensate collective mode fre- 
quencies and the density profile of the condensate at zero temperature |§] . To understand 
the nonzero temperature behaviour of Bose-condensed gases, several proposals have been 
made to generalize the Gross-Pitaevskii equation and to include the effects of the thermal 
cloud on the condensate. At the mean-field level this is achieved by introducing in the 
Gross-Pitaevskii equation real and imaginary terms, that describe the coherent and inco- 
herent effects of collisions between condensate and thermal atoms, and that in particular 
cause evaporation or growth of the condensate 0-0. However, at nonzero temperatures also 
fluctuations can play an important role. An example of this is the reversible formation of 
a condensate, as experimentally achieved by Stamper-Kurn et al. [|IIJ. Since the system is 
several times in the critical region, fluctuations of the order parameter around its mean-field 
value are of utmost importance to describe this experiment [n|. In addition, both quantum 



and thermal fluctuations are important to understand the stochastic nature of the collapse 
observed in 7 Li [J12 — 17}] and the phenomena of phase 'diffusion' [18|, in which case these 



fluctuations disturb the global phase of the condensate. Finally, from a fundamental point 
of view, a consistent description of a partially Bose-Einstein condensed gas requires that 
the fluctuation-dissipation theorem is obeyed, since this ensures relaxation of the system 
towards its correct physical equilibrium. Therefore, if dissipation is to be included in the 
generalized Gross-Pitaevskii equation, fluctuations must also be included. 

Gardiner and Zoller have included such fluctuations in the description of a Bose- 
condensed system by deriving with second-order perturbation theory a master equation 
for the one-body density matrix [ 19| , a procedure well-known from quantum optics. How- 



ever, in this paper we will use the nonperturbative formulation developed previously by 
one of us |f7|j20|]. Using field-theoretical techniques Stoof derived a Fokker-Planck equation 
describing the full nonequilibrium probability distribution of the order parameter. An equiv- 
alent formulation of this theory can be given in terms of a dissipative nonlinear Schrodinger 
equation with noise. Although in principle we can turn to numerical methods for the solu- 
tion of the Fokker-Planck equation, or its corresponding Langevin equation |TTJ , we find it 



more convenient here to proceed analytically, by means of a variational method. Variational 
approximations have previously provided a useful way to make analytical progress, and cap- 
ture as many of the physics as possible. In particular, when applied to the zero-temperature 
Gross-Pitaevskii equation, a gaussian variational approximation has led to good results on 



the collective modes of the condensate |21^p3|, and on the description of the macroscopic 



tunneling of a condensate with attractive interactions |11^J15|J24|J25|| . It is the aim of this 
paper to apply a similar variational method also to the dissipative nonlinear Schrodinger 
equation with noise appropriate for nonzero temperatures. We achieve this by assuming 



that the thermal cloud is in equilibrium at all times, and therefore acts as a 'heat bath' on 
the condensate. The stochastic nonlinear Schrodinger equation then obeys an equilibrium 
version of the fluctuation-dissipation theorem, which ensures that the condensate relaxes 
to the physically correct equilibrium. With this assumption, we are then able to derive 
Langevin equations for the variational parameters in our gaussian ansatz, which turn out 
to be equivalent to the equations of motion for a Brownian particle in a potential. These 
equations are coupled to a stochastic rate equation for the number of atoms in the conden- 
sate. Using these equations of motion, we are then able to describe collisional damping of 
the condensate collective modes at nonzero temperatures, and the condensate growth and 
stochastic initiation of the collapse, as recently observed in 7 Li [|T6|,|T7 . 



The rest of this paper is organized as follows. To make this paper selfcontained, we 
review in Sec. II the techniques of path integrals, and their application to stochastic dif- 
ferential equations. We use the method of functional integration throughout this paper. In 
Sec. Ill we review the Fokker-Planck equation describing the nonequilibrium dynamics of 
a Bose-Einstein condensed gas, and discuss the equilibrium solution of this Fokker-Planck 
equation. The most important result of this section is the Langevin field equation for the or- 
der parameter, that obeys the fluctuation-dissipation theorem. This Langevin field equation 
takes the form of a dissipative nonlinear Schrodinger equation with noise. We also derive the 
stochastic equations of motion for the density and phase of the condensate, and a damped 
wave equation describing the propagation of sound waves in a homogeneous Bose gas, at 
nonzero temperatures. In Sec. IV we present the variational approximation to our nonlinear 
dissipative Schodinger equation with noise, and also derive stochastic equations of motion 
for the variational parameters, which are the central result of this paper. To the best of our 
knowledge, a variational method for stochastic field equations, such as the stochastic non- 
linear Schrodinger equation under consideration here, has not been developed previously. 
In Sec. V we apply our equations to calculate the temperature dependence of the damping 
and frequencies of the collective modes of a condensate, and to obtain a description of a 
growth-collapse curve of a condensate with attractive interactions. We end in Sec. VI with 
our conclusions. 



II. PATH INTEGRALS AND STOCHASTIC DIFFERENTIAL EQUATIONS 

In this section, we discuss the application of path integrals to stochastic differential 
equations. First, we consider Brownian motion in the overdamped limit, which means that 
the acceleration of the particle can be neglected with respect its damping rate. In the second 
part we consider underdamped Brownian motion, in which the nonzero mass of the particle 
plays a crucial role. 



A. Overdamped Brownian motion 



Consider the one-dimensional Langevin equation fl26|,27l 

q(t)=f(Q(t))+v(t). (1) 

In this equation, q(t) denotes the position of the particle, and rj(t) is a gaussian noise term, 
which has a time- correlation function given by 



( V {t')r ] {t)) = c5(t'-t). (2) 

Here, the brackets denote averaging over different realizations of the noise. The parameter a 
is positive and real, and is often referred to as the 'strength' of the noise. For a given initial 
condition q(t ) = q , the Langevin equation in Eq. (P generates a probability distribution 
P[q,t;qo,t Q ]. By definition, P[q,t;q ,t ]dq denotes the probability that a solution q(t) of 
Eq. (0) starting at q at time t , has a value between q and q + dq at time t. From this 
definition, we can immediately write down an expression for this probability distribution, 



P[q,t;q ,t ] = (5(q(t)-q)). (3) 

Here q(t) denotes again a solution of the Langevin equation for a particular realization of 
the noise, starting at q(t ) = q . The brackets denote averaging over different realizations 
of the noise, as in Eq. (0). 



We now want to derive a path-integral expression |28| , [29| for the probability distribution 
P[q, t; qo, to]. To achieve this, we first divide the time interval t — t into N pieces, each of 
length A = (t — t )/N. Using the notation q(t n ) = q n , and rj(t n ) = rj n , we discretize the 
Langevin equation in Eq. ([I]) as follows 

^(?n+l - Qn) = f(q n ) + Tin- (4) 

The time-correlation of the noise is given by 

(ViVj) = ^<%, (5) 

which reduces to Eq. (g) in the limit A — > 0. Making use of the fact that the noise has a 
gaussian distribution, we now write the probability distribution for a particular realization 
of the noise as 



*({%» = (-5-) -p --E;^ ■ (6) 



Using this probability distribution, we are able to calculate averages over the noise as gaus- 
sian integrals. The two-point function of the noise is, for example, given by 



.N-l 
(ViVj) = / II d VnViVjP({Vn}), (7) 



which reproduces the correct correlations, given in Eq. @. We now first calculate the 
probability distribution P[qi,ti,qo,t Q ], which is the probability distribution that a solution 
of the Langevin equation in Eq. (HD reaches the value q\ at time t\ = to + A. Since we can 
easily solve the discrete Langevin equation in Eq. (^) explicitly for one time step, we find 
from the definition in Eq. (|3|) that 

P[qi,t 1 ;qo,t ]= J dr] (-£-) exp |-— rft \5 (q + A(/(g ) + Vo) ~ Qi) ■ (8) 



Secondly, we integrate out the noise r] to obtain the result 

P[quh;qo,t ] = ^xp j-A (Sl^. _ /(gb )) 2 J . (9) 

We then use the this expression at each time step, and 'tie' them together using 

P[ft+i,£i+i;g»_i,*i-i] = / dqiP[qi + i,ti +1 ;qi,ti]P[qi, **; gi_i, *»-i], (10) 

which follows from the fact that the total probability is conserved. The result for 
P[qN,tN',qo,to] then becomes, after a combination of Eqs. (U) and ( |10|) at each interme- 
diate time step, 

fN-l \ ( a JV-1 



Pfev, tiv; go, to] = A-(^ 2 ) / n <*& exp - g f; (^-^ - /(ft)) • (11) 



Note that this expression explicitly shows that the integration is only over intermediate 
coordinates, and that the boundary values q^ and go are fixed. We now take the limit 
N — > 00 and A — > 0, while keeping t]y — t fixed. If we absorb the prefactor in Eq. (]TT|) in 
the integral measure, we get, after putting q N = q and ijv = t, the result 

P[q,t;q ,t ]= d[q]e iS ^' n . (12) 

Jq(t )=q 

In the exponent of the integrant we extracted a factor i/h, where h is Planck's constant, 
for reasons that become clear shortly. The integral measure of the functional integral in 
Eq. ( |12|) now denotes integration over all paths q(t) with boundary conditions q(t ) = q Q , and 
q(t) = q. Each of these paths gives a weighted contribution to the probability distribution, 
with a weight factor proportional to e lS ^ q ^ h . The action S[q] is given by 

S [ q ] = l nfdt'^{ q {t')-f{ q {t'))f. (13) 

From the path-integral expression in Eq. ( |T^ ) we can now easily derive a partial differential 
equation for P[q, t; qo,to]- This equation will turn out to be the Fokker-Planck equation 
2T||2"T|] corresponding to the Langevin equation of Eq. ([!]). To derive the Fokker-Planck 
equation, we make a connection with quantum mechanics. 

From ordinary quantum mechanics for a point particle with mass m in a potential V(q), 
we know that a matrix element of the evolution operator can be represented as a path 
integral [p8 , p9| , similar to Eq. ([12]) . Denoting that matrix element by W(q,t;q ,t ) = 
(q\e~ l( - t ~ t °' ,H ^'^' n \qo) , where H(p,q) = p 2 /2m + V(q) is the Hamilton operator, this path 
integral is given by 

W(q,t-q ,t )= d[q)e iSC W n , (14) 

Jq(t )=qo 

with S cl [q] the classical action for the particle 

S cl [q] = fdti (^nq 2 {t>) - V{q{t'))) . (15) 
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On the other hand, the evolution operator also obeys the time-dependent Schrodinger equa- 
tion 

d 
ih—W(q,t;q ,t ) = H(p q ,q)W(q,t;q ,t ). (16) 

Here, H(p q ,q) is the hamiltonian corresponding to the classical action in Eq. (|15D , in the 
position representation. It is essential that H(p q , q) is normal ordered, i.e., that the momen- 
tum operators are positioned left to the position operators, for the path-integral expression 
in Eq. fll~4]) to be valid. We can now derive the Fokker-Planck equation by noting that 
Eq. flT2| ) can be interpreted as the path-integral expression for the time evolution operator 
with a 'classical' action given by 

S[q}= fdt'Lit'), (17) 

J to 

where L(t) is by definition the lagrangian. The desired Fokker-Planck equation is, therefore, 
exactly the Schrodinger equation corresponding to the action S[q]. To derive the latter, 
we have to quantize this action, which can be done in a straightforward manner. The 
momentum conjugate to the position variable q is given by 

P« s § = 7(* -/(*))■ ( 18 ) 

The hamiltonian thus becomes 

H(p q , q) = p q q ~ L(q, q) = —p\ + p q f(q). (19) 

Note that we have positioned the momentum operators left to the position operators, in 



order to obtain the correct correspondence with Eq. (12). Next, we impose canonical com- 
mutation relations [q,p q ] = i%, so we write in the position representation p q = —i%d/dq. 
The 'Schrodinger equation' for P[q,t;qo,to] is then given by 

m dP[q,t,q ,t ] = H{pqA) p [qA%Mi (20) 



which can be explicitly rewritten as 

dP[q,t] d 
dt dq 



l|- /(?) 



P[q,t). (21) 



This is indeed the well-known Fokker-Planck equation corresponding to the Langevin equa- 
tion in Eq. ([!]). It has a 'streaming' term linear in the derivative with respect to q, which 
represents the reversible part of the Langevin equation. The 'diffusion' term, quadratic in 
the derivatives, represents the irreversible stochastic behaviour. All Fokker-Planck equations 
that correspond to a Langevin equation with gaussian noise have this general structure. Note 
that in writing down the Fokker-Planck equation, we dropped the dependence of the prob- 
ability distribution on the initial conditions. This is because the form of the Fokker-Planck 
equation is independent of these initial conditions. Although the hamiltonian in Eq. ( |T§D 
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is not hermitian, the Fokker- Planck equation in Eq. ( pl|) conserves the normalization of the 
probability distribution, as can be seen from 



d_ 
It 



dqP[q,t] = dq 



d_ 
dq 



a_d_ 
2dq 



m 



P[q,t] 



a_d_ 
Ydq 



m 



P[q,t] 



0. 



(22) 



where we made use of the fact that a normalized probability distribution has to vanish at 
the boundaries. The equilibrium solution of the Fokker-Planck equation is given by 



P[q,t -»• oo] oc expi-y dq'f(q') 



(23) 



Since we have used units such that the mass of the Brownian particle and the friction 
coefficient are equal to one in the Langevin equation, we see that, in order for the probability 
distribution to relax to the correct Boltzmann distribution, we have to take the strength of 
the noise equal to o = 2k^I ', where k-Q is Boltzmann's constant, and T is the temperature. 
This is the fluctuation-dissipation theorem, which causes the probability distribution to relax 
to the correct physical equilibrium. 

At this point we would also like to make clear that in writing down the time sliced 
version of the Langevin equation, we have made the choice to interpret the noise term as a 
so-called Ito process, as opposed to a Stratonovich process. The difference between Ito and 
Stratonovich calculus emerges when one deals with multiplicative noise. For example, let us 
consider the Langevin equation 



q(t)=f(q(t))+g(q(t))r ] (t), 

where r)(t) is a gaussian noise term with correlations given by Eq. 
noise in Eq. (p4|) as an Ito process, the discretization reads 



(24) 
If one interprets the 



A 



(q n +i - q n ) = f(q n ) + g(q n )vn- 



The corresponding Fokker-Planck equation then becomes 



dP l [q,t] d 
dt dq 



Ufa - m 



p%t}. 



(25) 



(26) 



However, the time sliced version of Eq. fl2~j|) is in the Stratonovich calculus given by 



A 



(q n +i ~ q n ) = g [/(?»») + f(q n +i)] + g [0(9*0 + 9(Qn+i)]Vr> 



In this case, the Fokker-Planck equation reads 

0P s [q,t] 



dt 



d_ 

dq 



rr 2' - d 



2^- /(*) 



P%t}. 



(27) 



(28) 



The Stratonovich interpretation therefore leads to an additional noise-induced drift term in 
the equation for the average of q(t), as can be seen from 



^T^ = (f(<l)) S (t) + cr(g( q )g'(q)) s (t), (29) 

where g'(q) = dg/dq. This result follows straightforward from the Fokker-Planck equation 
in Eq. (|28|), with the use of partial integration. Note that the second term on the right-hand 
side of the last equation, which is the so-called spurious or noise-induced drift term, is absent 
in the case of an Ito process. In physics, a Stratonovich process arises naturally when the 
delta function in the time correlation of the noise is the result of a limiting procedure in 
which the correlation time becomes equal to zero. 

With these important remarks we conclude the discussion of the one-dimensional 
Langevin equation and its corresponding Fokker-Planck equation. As a further illustra- 
tion of the path-integral methods, we now discuss the case of the Brownian motion of a 
massive particle in a potential. 



B. Underdamped Brownian motion 

Let us consider the Langevin equation for the Brownian motion in one dimension of a 
particle with mass m in a potential V(q), 

q(t)+iq(t) = -^(q(t))+v(t). (30) 

In this equation 77 (£) is a fluctuating force per unit mass, with a gaussian probability dis- 
tribution. The parameter 7 > is a friction constant. The time-correlation of the noise is 
given by 

(v(t'Mt)) = ^(t'-t). (31) 

Here j3 = 1/ksT is the inverse temperature. Note that the strength of the fluctuations is 
related to the amount of dissipation 7 through Eq. fl3~T|). As mentioned previously, this is 
the fluctuation-dissipation theorem, which ensures that the probability distribution for the 
position and velocity of the Brownian particle relaxes to the Boltzmann distribution, as we 
will see in detail lateron. We now want to derive the Fokker-Planck equation associated 
with this Langevin equation of motion. To do so, we first write it as a set of two first-order 
differential equations 

q(t) = v(t); 

v(t) = -yv(t) - -^~(q(t)) + V(t). (32) 

m oq 

We are interested in the probability distribution P[q,v,t;q ,v ,t ], which is defined as the 
probability density that a particle with velocity Vq and position q at an initial time to, has 
a velocity v and a position q at time t. This probability distribution is thus given by 

P[q,v,t;q ,v ,t } = (6{q(t) - q)S(v(t) - v)), (33) 



where (q(t),v(t)) is a solution of Eq. (|3^), with initial conditions (q(to),v(to)) = (q ,vo). 
Following the same time-slicing procedure as in the previous section, we can derive a path- 
integral expression for this probability distribution. This path integral is in first instance 
given by 

P[q, v, t; g , v , t ] = / d[q] / d[v] S[v(t') - q(t')] 

^{-tf^{ iV) + lv(i) + ^ m) )\ (34) 

We next represent the delta functional by a Fourier path integral over an auxiliary coordinate 
p q . As the notation suggests, this turns out to be the momentum conjugate to q. The path- 
integral expression for P[q,v,t;qo,vo, to] then becomes 

P[q,v,t;q ,v ,t } = / d[q] / d[p q ] / d[v] 

Jq(t )=q J Jv(t )=v 

x exp U jTdf L{t'){q{t') - v(t')) + ^ ^(0 + 1V {t') + ^(Q(t'))J) | • (35) 

In this expression, we again extracted a factor of i/h outside the time integral in the expo- 
nent to make the connection with quantum mechanics more obvious, and we also absorbed 
some normalization factors in the path-integral measure. We then introduce the momentum 
conjugate to v, denoted by p v , by multiplying the integrant in Eq. (|35| ) by a factor one, 
written as the gaussian functional integral over p v given by 

1 =/d W expU/V^ U) - ^ (*(0 +7.(0 + ^WO))) 2 } ■ (36) 



This procedure is generally known as a Hubbard-Stratonovich transformation ||31|| . After 
this procedure, the result for the probability distribution P[q, v,t; qo,Vo,to] reads 

/•g( i )=9 f rv(t)=V r 

P[q,v,t;q ,vo,t ] = / d[q] / d[p q ] / d[v] / d[p v ] 

■ J q{to)=qo J Jv(t )=v J 

x exp |^ 4 di! (p q (t')q(t f ) + p v (t')v(t') - H{p q , q; p v , v))\ , (37) 

with a hamiltonian given by 

H(p q ,q; Pv ,v)= Pq v - ^ - p v (t« + ~) • (38) 

At this point it might be somewhat confusing that the momentum conjugate to q is not 
simply proportional to v. We are, however, not quantizing a classical system, but instead 
trying to derive a path integral for the probability distribution generated by a classical 
stochastic equation of motion. The connection with quantum mechanics lies in the fact that 
we can identify this probability distribution with a quantum mechanical amplitude for some 
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quantum system. This does not mean that the Brownian particle has wave-like properties. 
Nevertheless, Eq. (^) is precisely the canonical path-integral representation for a matrix 
element of the evolution operator. Therefore, the probability distribution P[q,v,t;q ,vo,t ] 
obeys the time-dependent Schrodinger equation with the hamiltonian H(p q ,q;p v ,v) given 
by Eq. fl38|), in the position representation. It is again essential that we use normal ordering. 
We can thus quantize this hamiltonian, by putting [q,p q ] = [v,p v ] = i% with all other 
commutators equal to zero. So we have in the position representation p q = —ihd/dq and 
p v = —ihd/dv. The Schrodinger equation now results in 



dP[q,v,t] 
dt 



d d ( 1 dV\ 7 d 2 



dq dv \ m dq J m/3 dv" 1 



P[q,v,t], (39) 



which is indeed the Fokker-Planck equation associated with the Brownian motion of a parti- 



cle with mass m in a potential V(q), and known as the Kramers-Klein equation fl32|. It can 



be shown that the stationary solution of Eq. (|39|) is given by the Boltzmann distribution 

P[q, v,t^oo]oc exp 1-/3 (-mv 2 + V(q)\ } , (40) 

which may be checked by insertion. It is important to realize that the fluctuation-dissipation 
theorem in Eq. (|3T| ) is again essential for the probability distribution to relax to the correct 
equilibrium distribution. It embodies the fact that dissipation and thermal fluctuations 
cooperate to achieve thermal equilibrium. This concludes our brief review of the connection 
between functional integration and stochastic differential equations. In the next section we 
use these techniques in the treatment of the nonequilibrium dynamics of a Bose-Einstein 
condensate. 

III. NONEQUILIBRIUM DYNAMICS 

In this section, we present the Fokker-Planck equation describing the nonequilibrium 
dynamics of a Bose-Einstein condensed gas, and its corresponding Langevin field equation. 
The so-called hydrodynamic formulation will also be discussed. Since the Langevin field 
equation generalizes the Gross-Pitaevskii equation to nonzero temperature, we start our 
discussion by recalling this well-known equation. 

A. Stochastic nonlinear Schrodinger equation 

The dynamics of a trapped Bose-Einstein condensate is at sufficiently low temperatures 
very well described by the time-dependent Gross-Pitaevskii equation 

ih ^ajr = {-^r + ^ cxt (x) +T 2B i*(x,t)i 2 | *( X ,t), (4i) 

where h is Planck's constant, m is the mass of a single atom, V ext (x) is the external trapping 
potential and T 2B = Airah 2 /m is the two-body transition matrix, with a the s-wave scattering 
length. The Gross-Pitaevskii equation arises as the equation of motion for the superfluid 
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order parameter, which is the expectation value of the Bose field operator ^(x, t), that 
annihilates an atom at position x and at time t. The Gross-Pitaevskii equation is also 
referred to as the non-linear Schrodinger equation for the macroscopic condensate wave 
function \l/(x, t), since the condensate density is given by 

n c (x,£) = |^(x,£)| 2 . (42) 

The time- dependent Gross-Pitaevskii equation has stationary solutions of the form \I/(x, t) = 
, I'(x)e JA '' //?i , where the parameter /i is the chemical potential that fixes the number of atoms 
in the condensate and ^(x) now obeys the time-independent Gross-Pitaevskii equation, 

1 " V2 +y cxt (x)-/i + T 2B |$(x)| 2 W(x) = 0. (43) 



[ 2m 

The time- dependent Gross-Pitaevskii equation is a semiclassical mean-field equation, 
describing the average dynamics of the condensate only. It contains no description of the 
relaxation of the condensate towards equilibrium, and neither does it contain condensate 
growth from the thermal cloud, or condensate evaporation, at nonzero temperatures. More- 
over, it completely neglects fluctuations of the order parameter around its mean value in 
the description. Therefore, we would like to modify the Gross-Pitaevskii equation such that 
it contains fluctuations due to incoherent collisions between condensate and noncondensate 
atoms, as well as condensate growth and evaporation. In order to do so consistently, we 
have to consider the full probability distribution for the order parameter, which can be found 
by means of the many-body T-matrix approximation to a field-theoretic formulation of the 
Keldysh theory |7,20]. It is given as a functional integral by 



P[<f>, 0*; t] = [* ^ * W d[<f>*]d[<j>] exp \ % -S^ [0*, 0]) , (44) 



with an effective action 
S eff [0*,0] = f dt' fdx 

J t 



ft£ K (x,£') 



X 



' %n h + lS~ ~ ^ ext(x) + ^ x > t ') + MO + T 2B |0(x,t')l 2 ) 0( x ,t'; 



• (45) 



In this effective action, the imaginary term iR(x, t) describes the exchange of atoms between 
the condensate and the thermal cloud. Since, at this point, we also want to be able to 
describe a thermal cloud that is not in thermal equilibrium, we have to allow for a time 
dependent chemical potential. Before we discuss the physical content of the expressions 
in Eqs. (|44] ) and (|45| ) further, let us first derive the Fokker-Planck equation determining 
the time-dependence of P[<p,<fi*;i\. To do so, we note that the expressions in Eqs. (f44"|) 
and ( f4"5|) are very similar to the path-integral expressions we encountered in the previous 
section for the probability distribution generated by a stochastic differential equation. The 
only difference between Eqs. QT2] ) and ([P]) is that the functional integration is now over 
all complex fields 0*(x, t) and 0(x, t), instead of real functions q(t). Also note that we did 
not specify the initial conditions at the time to- This is because we are only interested in 
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the universal long-time dynamics of the gas, which are independent of the specific form of 
the initial conditions. Moreover, as we have seen in the previous section, the form of the 
Fokker-Planck equation is in fact independent of these initial conditions. From Eq. ( f44|) we 



can derive the Fokker-Planck equation by quantizing the effective action in Eq. ([45]) , just as 
in the previous section. It is ultimately given by 




h 2 v 2 

2m 

n 2 v 2 



- V oxt (x) - n{t) - zi?(x, t) + T 2B 
+ V ext {x)-v{t)+tR{x,t) + T 2B 




2m 

fr£ K (x,£)P[0*, </>;£]. 

This Fokker-Planck equation describes the time evolution of the probability distribution of 
the condensate wave function at nonzero temperatures, in the presence of a thermal cloud. 
The dissipation term i?(x, t) describes the exchange of atoms between the thermal cloud 
and the condensate, due to elastic collisions. In the Hartree-Fock approximation, which is 
sufficiently accurate for the nonzero temperatures of interest here, it is given by p0| 



fl ,x,^ 2 ,(T»)7^L/*/^ (2 , mkl _ k2 _ k3 



(2tt) 3 J (2vr) 3 J (2tt) 
xd (e + ei - e 2 - e 3 ) [^(1 + N 2 )(l + N 3 ) - (1 + NJN^} , (47) 

In this expression, N{ = iV(x, kj, t) is the Wigner distribution function of the thermal cloud, 
which can be determined by solving the corresponding quantum Boltzmann equation. We 
will not do this explicitly here, since lateron we assume that the noncondensed cloud is in 
thermal equilibrium. The energy of a thermal atom is given by 

6i = ^if + yCXt(x) + 2T2B I ^WX*)! 2 - ( 48 ) 

Note that both in this expression, and in the Fokker-Planck equation in Eq. ( |46f) we neglected 
the effect of the mean field of the thermal atoms, because it plays a minor role in the dynamics 
of the condensate. Note also that, for the average value of the order parameter calculated 
with the probability distribution in Eq. ((44]), we used the notation (</>(x))(£). The noisy 
order parameter field will be denoted by </>(x, t) and for stochastic averages of this quantity 
we will use the notation (</>(x, £)). Since the Fokker-Planck equation and its corresponding 
Langevin equation are equivalent, we have of course that (</>(x))(£) = (</>(x, t)). The Keldysh 
self-energy fr£ K (x, t) in the Fokker-Planck equation describes the thermal fluctuations due 
to incoherent collisions between condensate and noncondensate atoms. It is given explicitly 



by [20 



^M-«n7^/^/|?w* 



x5 (e + ei - e 2 - e 3 ) [JVi(l + N 2 )(l + N 3 ) + (1 + N 1 )N 2 N 3 ] . (49) 
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Note that both the dissipation i?(x, t) and the Keldysh self-energy frE K (x, t), depend on the 
energy e to take a condensate atom out of the gas at position x and time t, which has to be 
determined selfconsistently. This implies that e is actually an operator in the configuration 
space of the order parameter, and given by |2(| 



ft 2 V 2 

2m 



+ 1/ ext (x)+T 



2B| 



(50) 



The fact that e should be viewed as an operator will turn out to be crucial for the probability 
distribution of the order parameter to relax to the correct equilibrium distribution function. 
Although our Fokker-Planck equation for the condensate, coupled to the appropriate 
quantum Boltzmann equation for the Wigner distribution function of the thermal cloud, 
describes in principle the full nonequilibrium dynamics of the Bose-condensed gas, its solu- 
tion is very difficult even numerically. This is because of the fact that in the Fokker-Planck 
equation the dissipation -R(x, t) and the Keldysh self-energy also depend on the condensate 
wave function, through their dependence on e given in Eq. (|50|), and through the mean- field 
effect of the condensate on the thermal atoms. As a result, writing down the correspond- 
ing Langevin equation results in a stochastic equation with multiplicative noise, and with 
a prefactor of the noise that has a complicated dependence on 0(x, t). We can, however, 
make progress by assuming that the thermal cloud is sufficiently close to equilibrium, which 
is, for example, justified for linear-response calculations around equilibrium, and also for 
condensate growth if the evaporative cooling is performed sufficiently slowly. From now on, 
we therefore assume that the thermal cloud can be described by a Bose distribution function 



N& 



,/3(«i-A*) 



1 



-1 



(51) 



with a chemical potential /i and an inverse temperature f3 = l/k^T. The thermal cloud 
therefore now acts as a 'heat bath' on the condensate. Making the e-dependence explicit for 
a moment, we can relate the dissipation i?(x; e), and the Keldysh self-energy frE K (x; e) by 
means of 



zi?(x;e) = -^E K (x;e)[H-2AT(e)]- 1 , 



(52) 



which follows simply from the form of the Bose distribution function, together with the 
energy conserving delta function in Eqs. (|47|) and fl49f ). This relation between the dissipation 
i?(x; e) and the Keldysh self-energy fi£ K (x; e) determining the strength of the fluctuations, 
is in fact the fluctuation-dissipation theorem. Just as in the case of the Brownian motion 
of a particle discussed in the previous section, it causes the system to relax to the correct 
equilibrium distribution, as we will see below. Since we are dealing with Bose condensation, 
the occupation numbers N(e) are generally very large, and we have in a good approximation 



[l + 2N(e)Y 



A*)]- 



(53) 



If we combine this result with Eq. 
at the approximation 



2|), and substitute the operator in Eq. 



we arrive 



iR(x,t) ^-jhE 



K, 



X,t) 



h 2 v 7 

2m 



V 



cxt / 



fi + T 



2B| 



(54) 



13 



where fr£ K (x, t) = fr£ K (x; (/U c (x, £))), and the local chemical potential of the condensate 
<> c (x,£)) is given by 

WM)> = w(x))(t)|2 /^x(0*(x))(t) (-^21 + V cxt (x) + — K0(x)>(i)| 2 J (0(x))(t). 

(55) 

We now show that the above 'classical' approximation to the fluctuation-dissipation 
theorem indeed leads to the correct equilibrium. Let us therefore substitute Eq. (0) into 
the Fokker-Planck equation, which simplifies to 

-| / dx ft£ K (x, t)j^ ("^ + V cxt (x) -im + T 2B |0(x)| 2 ^) 0(x)P[«f , 0; t] 
-£ / rfx ft£ K (x, t)-^ (-^ + V cxt (x) -ii + T 2B |0(x)| 2 ) 0*(x)P[0*, 0; t] 

-5/ & ^l?w' aiK(x '* W ^ fl - (56) 

The stationary solution of this Fokker-Planck equation is given by 

,2v72 ^2B \ 1 



P[<f>", </>; t -► oo] oc exp 1-/3 1 rfx 0*(x) f -^- + ^ cxt (x) - fi + 



x)| 2 U(x) , (57) 



as can be checked by substitution. To see that Eq. fl5"T|) is in fact the correct equilibrium 
distribution, we have to show that the macroscopic condensate wave function (0(x)) obeys 
the time-independent Gross-Pitaevskii equation. To see this, we first note that 

/d[0*]c#]^^P[0*,0;£]=O (58) 

for a general probability distribution which vanishes at the boundaries of the domain of 
integration. If we apply this to the equilibrium distribution P[|0|;£ — ► oo] we get, by 
applying the mean- field approximation (|0(x)| 2 0(x)) ~ |(0(x))| 2 (0(x)), the desired result 

' /r ^+\/ ext (x)-/i + T 2B |(0(x))| 2 V0(x)) = O, (59) 



Ira 



which is precisely the time- independent Gross-Pitaevskii equation. Note that Eq. (54) 
together with the time-independent Gross-Pitaevskii equation implies that in equilibrium 
(P(x, t)) = 0. This means that there is detailed balance between the condensate and the 
thermal cloud, and that there is, on average, no condensate growth or evaporation when the 
system has relaxed to equilibrium. 

Using the results of the previous section, we now give a formulation of the nonequilibrium 
theory discussed above in terms of a Langevin field equation corresponding to the Fokker- 
Planck equation in Eq. (0). This Langevin field equation takes the form of a dissipative 
nonlinear Schrodinger equation with noise, given by 
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ih 



dct>(y:,t) 
dt 




h 2 v 7 

2m 



+ V cxt (x)- / i + T^|0(x,t)|^ U(x,t)+77(x,t). 



(60) 



This Langevin equation quite generally generalizes the Gross-Pitaevskii equation to nonzero 
temperatures, and includes both dissipation and thermal fluctuations. The complex gaussian 



noise in the Langevin field equation has correlations |20| . |TT 

ih 2 



(77*(x,t)77(x',t0) 



-E K (x,t)5(t-t / )^(x-x / 



(61) 



where the strength of the noise is determined by a Keldysh self-energy, given by 
SS k (x , „ = -« (T-)7 * / * / * (*)., (k, - k2 - k3 ) 



(2tt) 3 
x<5((/i c (x,t)) +ei -e 2 



(2vr) 3 y (2tt) 3 
e 3 ) [iVi(l + N 2 )(l + N 3 ) + (1 + iV^iVsiVg] 



(62) 



In this expression, iVj is again the Bose-distribution function of the thermal cloud, evaluated 
at an energy of a thermal particle, which is in the Hartree-Fock approximation given by 



2m 



+ V 



cxt / 



2T™\(<f>(x,t))\- 



(63) 



The average local chemical potential of the condensate atoms (/x c (x, t)), is given by Eq. (p5|). 
Notice that in the latter equation, and in the above expression for the energy of a thermal 
particle, (0(x, £)) has to be determined self-consistently, since only then the probability dis- 
tribution generated by the Langevin equation in Eq. fl60|) relaxes to the correct equilibrium. 
The stochastic non-linear Schrodinger equation in Eq. (|60"D , together with the expression 
for the Keldysh self-energy in Eq. fl6"2"|), gives a nonequilibrium description of the condensate, 
that obeys the fluctuation-dissipation theorem. In Sec. IV we use a variational ansatz 
to solve this equation. However, we first derive the corresponding noisy hydrodynamic 
formulation. 



B. Stochastic hydrodynamics 

The condensate is often described in terms of its density and its phase, by making the 
transformation <fi = ^fpe %e . When applied to the Gross-Pitaevskii equation, this transforma- 
tion results in the so-called Josephson equation for the phase, and a continuity equation for 
the density ||||. We now want to derive the generalization of these two equations to the 
case of our Langevin equation for the condensate. To do this, we first substitute the 'classi- 
cal' approximation to the fluctuation-dissipation theorem into the effective action S [<fr*, 4>], 
which now reads 



S' 



effr i * 



t f 2 

dt' / o?x — — 

t J ft£ K (x,£' 



2V72 



K 2 V 
2m 



V 



cxt 




(64) 



15 



The reason for this substitution is that we have defined the fluctuation-dissipation theorem 
as an operator equation in the configuration space of the order parameter, and not in terms 
of its density and phase. We can now easily substitute = ^/pe %e into this effective action. 
This substitution results in an effective action in terms of the density and the phase of the 
order parameter, i.e., 

S cS [p,6] = f dt' f <hc 

Jtn J 



ft£ K (x,t' 



\ dt' 4 v ' ; 2p(x,t') 



2 



+// c (x,t') -fi 



^P^ + v-W^O) 



+ ^S K (x,t')(Mc(x,t')-^P(x,t')') }, 



(65) 



where we used that S K (x, i) only has a negative imaginary part, as seen from Eq. flB"2|), and 
thus iS K (x, t) is a positive and real quantity. We defined the superfluid velocity v s (x, t), 
and the condensate chemical potential /i c (x, t) by means of Hffil , 



/x c (x, t) = J V + V cxt (x) + T 2B p(x, t) + -mv s 2 (x, t); 



2m^/p(x,t) 

v s (x,t) = -V0(x,t), (66) 

m 

which coincides with the expression in Eq. ([55]). The effective action S 1 [p, #] yields two 
stochastic equations of motion. The equation for the phase of the condensate takes the form 
of a stochastic Josephson equation 

t dfl(x,t) 13. ^ K( ^ h 2 V ■ (p(x,t)Vfl(x,t)) z/(x,t) 

*-fc 4^ E (x ' t} 2^7) = ^ _/ic(X ' ) + 7p^)' (67) 

Here, the real gaussian noise z/(x, t) has correlations given by 

ih 2 
(i/(x, t)z/(x', £')> = ^-fr£ K (x, *)<*(* - *')<*(* - x'). (68) 

The stochastic Josephson equation has two modifications with respect to the ordinary 
Josephson equation. First, it has a spatial diffusion-like term proportional to i^E K (x, t). 
This term will cause the phase to undergo spatial diffusion due to collisions of thermal 
atoms with the condensate atoms, not to be confused with the phenomenon of phase 'dif- 
fusion', which corresponds to spreading of the global phase due to quantum fluctuations 
18] , and therefore relax to a state where the phase is position independent. So, in equilib- 



rium we have (v s ) = 7i(V0)/m = 0, as expected. We will see lateron that this tendency 
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towards equilibrium will give rise to an increase in the sound velocity in the Bose conden- 
sate. Secondly, the Josephson equation has a noise term inversely proportional to the square 
root of the density. This noise represents the fluctuations in the phase of the condensate 
due to incoherent collisions of thermal atoms with condensate atoms, i.e., due to thermal 
fluctuations. 

The equation of motion for the density is a stochastic continuity equation with a source 
term, 

(J/i(X ' /) +V-(p(x,t)v s (x,t)) 



at 



-^S K (x, t) (p c (x, t) - p) p(x, t) + 2 v /p(x,£)£(x, t), (69) 



with correlations of the gaussian noise £(x, t) given by 

<£(x, t)£(x', t')} = ^M^(x - *)6(t - t>). (70) 

In appendix A it is explained that we have to interpret this noise as a Stratonovich process. 
This gives rise to additional drift terms in the equation of motion for the average of the 
density, because in Eq. ( p9| ) we are dealing with multiplicative noise. Note that from Eq. ( |69"D 
it is explicitly seen that there is condensate growth if \x > p c , i.e., if the chemical potential lies 
above the chemical potential of the condensate. If p < p c , there is condensate evaporation. 
We will omit here the Fokker-Planck equation in terms of p and 6, but only discuss 
the equilibrium distribution generated by Eqs. (^) and (|B"9"|). It is simply determined from 
P[(p*,(f);t — > oo] in Eq. fl57| ) by the substitution = ^fpe %i ', since the jacobian of this 
transformation is equal to one. So we have 



I f I h 2 V\ /p(x) T 2B 1 „ 

P[p, 9; t - oo] ex exp -(3 / rfxp(x) j^J- + V cxt (x) + — - p(x) + -mv s 2 (x) 

I J \ 2mJp(x) 2 2 



We see from this probability distribution that (v s ) = 0, which should be the case in equi- 
librium. The average density profile is again determined by the time-independent Gross- 
Pitaevskii equation, as explained before. 

To discuss the physical content of the stochastic continuity equation and the stochastic 
Josephson equation further, we now derive the wave equation describing the propagation of 
sound waves in a Bose Einstein condensate. For simplicity, we discuss here the homogeneous 
case, where \^ ext (x) = 0. A treatment of the trapped case is presented in Sec.V A. We 
linearize the averages of Eqs. ( |69"D and fl67|) around their equilibrium solutions (p(x, £)) = p 
and (v s (x, i)) = 0. Therefore, we write (p(x, £)) = p + 5p(x.,t), and (v s (x, £)) = 5v s (x, t), 
and substitute this into the average of Eqs. fl6"7|) and(|69|). Linearization results in two coupled 
equations of motion for the deviations, i.e., 

d6v 8 (x,t) 2B /3ih 2 T, K 2 

m — = -T V(<5p(x,t)) + — - — V (5v s (x,t)); 

^^+Po V ■ (5v s (x, t)) = -^S K (2T 2B p - p) 5p(x, t). (71) 
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Note that we have made use of the fact that frX K is independent of the spatial coordinates 
for a homogeneous Bose gas, as can be seen from Eq. (0). Next, we combine these two 
equations to obtain a single damped wave equation for the propagation of sound waves in a 
homogeneous Bose gas, 

- C VJ«M) = -'^. (72) 

The relaxation time r is defined to be inversely proportional to the damping rate of the waves. 
Physically, this damping arises because the excitation of a sound wave slightly disturbs the 
equilibrium situation where the average growth or evaporation of the condensate is equal 
to zero. Hence, there is no longer detailed balance between the condensate and the thermal 
cloud, and the collisions between the condensate and thermal atoms drive the condensate 
back to the equilibrium situation, where <5p(x, t) = 0, and 5v s (x, t) — 0. The relaxation time 
r is given by 

- = ^S K T 2B p , (73) 

T Z 

where we used the time-independent Gross-Pitaevskii equation which reduces in this case to 
/i = T 2B po, to eliminate the chemical potential. The sound velocity c in Eq. (ff^) is given by 



c 2 



l + ^(/^ 2 



(74) 



where Cq = (T 2B p /m) 1 ^ 2 is the well-known zero temperature sound velocity, predicted by the 
Gross-Pitaevskii equation, and first obtained by Bogoliubov [[34]]. We see that our nonequi- 
librium treatment results in an increased sound velocity. This increase is a result from the 
term in the stochastic Josephson equation in Eq. (|67|) proportional to iKL K . Physically, this 
term represents the fact that the phase of the condensate undergoes spatial diffusion due 
to collisions between condensate and thermal atoms, and therefore relaxes to a state where 
(v s ) = 0. The spatial diffusion of the phase therefore increases the 'stiffness' of the conden- 
sate, and hence results in an increase of the sound velocity. Since the increase in the sound 
velocity is of order (9(|/5/iS K | 2 ), its effect is in general small below the critical temperature, 
as the collisionless limit is determined by |/3fo£ K | <C 1 and experiments are usually in this 
limit. The damped wave equation in Eq. (^) should be compared to the result found by 



Williams and Griffin [33|]. These authors use a dissipative non-linear Schodinger equation, 
with a damping term similar to Eq. (f|7]), to arrive at a damped wave equation describing 
the propagation of sound in a trapped Bose-Einstein condensate in the presence of a static 
thermal cloud. Note, however, that although the microscopic expression used by Williams 
and Griffin is of the same form as in Eq. (pT7|), the chemical potential of the condensate used 
by these authors in the calculation of i?(x, t) is not the operator given by Eq. fl50|). Instead, 
they use for the condensate energy in the energy- conserving delta function the expression 
(p c (x, £)), which in principle violates the fluctuation-dissipation theorem. As a consequence, 
these authors do not find an increase in the sound velocity at nonzero temperatures. 



IV. VARIATIONAL APPROXIMATION 

Although the stochastic non-linear Schrodinger equation given in Eq. (|SUD, or equiva- 
lently, the hydrodynamic formulation given in Eqs. (0) and (|69|), give a full nonequilibrium 
description of the condensate that can in principle be solved numerically JO]] , we find it more 
convenient to make analytical progress. Therefore, in the case of a harmonic trapping po- 
tential V^ ext (x) = J2j muJ j x 'j/ 1 2, we consider a gaussian variational ansatz for the condensate 
wave function 

«x, t) = ^je*(" n y^j) ' exp {-^L (, _ "L^ (t) ) } . (75) 

Here, the variational parameters ?j(t) denote the gaussian widths of the condensate in the 
three spatial directions. The wave function is normalized to the number of atoms in the 
condensate N c (t). This ansatz is different from the ones used in previous work, in the sense 
that it also contains a global phase 9o(t). This turns out to be crucial, since the number 
of particles N c (t), which is the variable conjugate to the global phase, is not constant in 
our case. Therefore, one must also allow for fluctuations in the global phase 9 (t) of the 
condensate. We expect this ansatz to give correct results when the number of particles 
is small, because the mean-field interaction of the condensate will then be small, and the 
condensate density profile will be close to the ideal gas solution. Moreover, it has also proven 
to give correct results for the frequencies of the collective modes of the condensate, even in 
the Thomas-Fermi regime, where the mean-field interaction, and thus the number of atoms 



in the condensate, is large [pi]] . Therefore, we expect also to obtain physically sensible results 
even in this case. 

When the gaussian ansatz is applied to the Gross-Pitaevskii equation we find that the 
variational parameters qj(t) obey Newton's equations of motion |13|,2i,24 



1 dV 

-mN c (t)q)(t) = -gj-Ht), N c (t)), (76) 



with a potential energy equal to 



W'ES + >^t £_"_'__ . (77) 



N c Ti 2 1 Ar 2 2 \ ah 2 N 2 c 

—5-5 + -mN c cu 2 q 2 + —= 2- 

Amqj 4 J J ^2nmq x q y q z 

From the action in Eq. (0), we want to derive similar equations of motion, extended to 
the nonequilibrium case. For simplicity, we first consider an ideal gas, i.e., we drop the 
mean-field interaction term T 2B |0(x, t)\ 2 . The condensate remains, however, in contact with 
the thermal cloud, that acts as a 'heat bath'. Secondly, we assume that the Keldysh self- 
energy is constant over the size of the condensate. Although this assumption is not justified 
in general, we can always approximately compensate for this, by calculating a position 
independent Keldysh self-energy hH K (t) by means of an appropriately averaged fr£ K (x, t) 
over the size of the condensate. We substitute our trial wave function into the effective 
action in Eq. fl64|) , to obtain a probability distribution in terms of -/V c , #o and q. It is given 
by 
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P[N c ,6 ,q;t]=Jd[N c ] d[9 ] d[q]exp{^ cff [iV c A,q]}, (78) 

with an effective action that reads 

s cS [n c , 9 0] q ] = £ d^—jjv^) L^p- + ^(0 - E ^(0 + E i"H&(0«b(0 - ^ 

+ £jg)e« 

+0(((3ihZ K ) 2 )\. (79) 

The potential in this effective action is defined by 



^*>-E@H 



mN c uffi\ , (80) 

which is precisely the potential given by Eq. (|77|) , without the mean-field interaction term. 
The condensate chemical potential for the gaussian ansatz is given by 

*(*) = |£(q(t), tfcW) + E ^J(t), (si) 

as expected. We now assume the dimensionless parameter ftihH K to be small, and thus 
restrict ourselves to a temperature regime sufficiently far below the critical temperature, 
where |/3/i£ K | <C 1, i.e., the collisionless regime. We can then to a good approximation 
neglect the terms quadratic in [3iTiYl K . The effective action in Eq. (|79|) thus becomes a 
sum of three squares, and we can extract equations of motion with gaussian noise terms, 
exactly as in Sec. II. Since the action is quadratic in 9 (t), we can integrate over this global 
phase exactly, because it only requires a gaussian integral. However, before we perform this 
integration, we discuss the stochastic equation of motion for 9 Q (t). With the techniques 
discussed in the Sec. II, we easily see that it is given by 

n d9^)_ = ^_ + 1 _ 1 + K£) 

dt i 4 i 4 >JN c (t) 

with time correlation of the noise 

(v{t')v{t)) = l -^±5{t'-t). (83) 

This stochastic equation again has the form of a Josephson equation with a noise term 
added, similar to Eq. (0). Note that the noise term in Eq. fl8"2|) is inversely proportional to 
the square root of the number of particles in the condensate. As a result, the Fokker-Planck 
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equation for the probability distribution of &o(t), associated with the Langevin equation in 
Eq. ( p!q) , will have a 'diffusion' term inversely proportional to the number of particles. This 
means that the global phase is only well determined if there is a infinite number of atoms in 
the condensate, otherwise the global phase undergoes phase diffusion, due to thermal fluctu- 
ations. This mechanism for phase diffusion is different than the phase 'diffusion' considered 
by Lewenstein and You |0|], who considered phase spreading due to quantum fluctuations. 
Having made these remarks, we perform the integration over 6o(t), and are left with a 
probability distribution for N c and q. It is given by 

P[N c ,q, v;t] = Jd[N c ] d[q] d[v] 6[v(t) - q(t')] exp {^ eff [iV c , q, v]} . (84) 

Here we introduced the velocity v(t) = q(t) by means of a delta functional. The resulting 
effective action reads 

^•■*'i - t«idmtim par + f !EK(i '» w 




Eqs. (|8~4"D and (|85|) are similar to the path-integral expressions we encountered in our discus- 
sion of the Brownian motion of a particle in a potential in Sec. II. Therefore we immediately 
conclude that the equations of motion for the variational parameters are given by 



l,nN Mm + MtlmWlU = -g(q(«), *(«)) + =&«, (86) 

with the time correlations of the gaussian noise terms £,j(t) given by 

jh 2 Y K (t'\ 
UtW)) = jL 1 A1 5 3k S(t-f). (87) 

So we have found the important result that the variational parameters obey the equations 
of motion of a Brownian particle with mass mN c /2 in a potential V(q, N c ). Physically, the 
variational description of the condensate with the Langevin equation in Eq. (|55|) , as opposed 
to Eq. ([76|) , has two important extra features. First, there is a damping term present, i.e., 
a term proportional to the velocity (jj(t). This damping term can, for example, be used 
to calculate the damping on the collective modes of the condensate. Since the damping 
term is proportional to /i£ K , we conclude that it arises because of incoherent collisions 
between condensate and thermal atoms, which drive the condensate back to equilibrium, 
and let the phase of the condensate relax to a state where the phase is, on average, position 
independent. Secondly, since the Langevin equation also contains fluctuations, it can, for 
example, be used to describe the stochastic initiation of the collapse observed in 7 Li fl4HT7 . 
Our description contains thermal fluctuations, that cause the condensate to overcome the 
macroscopic energy barrier and start the collapse. In the next section, we will present the 
result of calculations that we have done on the two above mentioned phenomena. 

Since the potential in the Langevin equation in Eq. fl86|) depends on the number of 
condensate particles N c , we have to couple Eq. (|86|) to a rate equation for the number of 
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J5) 



condensate atoms. This stochastic rate equation also follows directly from the effective 
action in Eq. (jSq) with the techniques discussed previously, and is given by 



dN c (t) 
dt 



*£ K (£) \fjL c {t) - /i] N c {t) + 2^Ut)r](t), 



with correlations of the gaussian noise given by 



(v(t'Mt)) = ^^5(t'-t). 



(88) 



(89) 



As explained in appendix A, we have to treat the multiplicative noise in Eq. ( JS8|) again as 
a Stratonovich process, to achieve the correct equilibrium distribution. Physically, Eq. ( 58 ) 
describes the growth or evaporation of the condensate. The noise term in Eq. ( |S8D represents 
the fluctuations in the number of particles. 

To see that Eqs. (|86| ) and ( |88| ) generate the correct equilibrium distribution, we now 
discuss the Fokker-Planck equation for P[N C , q, v;t]. However, let us first discuss the equi- 
librium solution we expect on basis of Eq. (|57|). A substitution of the gaussian ansatz in 
Eq. ((75|) into equation Eq. (|57|) results in 



P[N C , q, v; t — > oo] oc exp < — /3 



1 



*E^mN c v] + V(q,N c ) - /iN c 



(90) 



where V(q,N c ) is the potential given by Eq. (|77|). Although the Langevin equations for 
qj(i) and N c (t) did, in first instance, not include the mean- field interactions, we argue that 
they also are correct for the interacting case. The reason for this is, that in this manner we 
are led to the correct equilibrium distribution as we show now. Let us therefore determine 
the Fokker-Planck equation for the probability distribution of N c , q and v, generated by 
the stochastic equations in Eqs. (J36D and (|88|) with the interacting potential in Eq. ([77]). It 
is given by 



dP[N c ,q,v;t] 



E 



dt 
d d fpih 2 E K (t) 



d 

dN r 



dv 



-zS K (t) 



j v 2mq'j 



2 dV 
mN c dqj 



[q,N c )) + 



i^S K (t) d 2 
m 2 N c q 2 dv 2 



dV 

dN r 



^ N c ) + J2 - A mv) - fi 



A , , tZ K (t) d Ar d 
N c + — _ v ; — N r 



2 dN r dN c 



P[iV c ,q,v;t] + 

P[iV c ,q,v;t], (91) 



and insertion of the equilibrium distribution shows that it is indeed a stationary solution 
of this Fokker-Planck equation. Thus we conclude that the Langevin equations for qj(t) 
and N c (t) give, with the potential V(q, AT C ) given in Eq. (f77[) , the correct description of 
the nonequilibrium dynamics of a Bose-Einstein condensate in the gaussian approximation. 
This description includes damping of the collective modes of the condensate, as well as 
condensate growth and evaporation. The essence of our method lies in the fluctuation- 
dissipation theorem, which ensures the relaxation towards the correct physical equilibrium 
distribution. 
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V. APPLICATIONS 

In this section we first apply the Langevin equations for the variational parameters, 
derived in the previous section, to the calculation of the damping and frequency of the 
collective modes of the condensate. As a second application, we also obtain a description of 
the initial growth of a condensate. 

A. Collective modes of the condensate 



In this section we use the Langevin equations in Eq. ( j86|) for the gaussian variational 
parameters, and the stochastic rate equation in Eq. ([88]) for the number of particles in the 
condensate, to obtain a description of the collective modes of the condensate. We calculate 
the frequency and damping of both the monopole and quadrupole mode in an isotropic trap, 
and compare with the theoretical results found by Williams and Griffin p5|p6| . Since we are 
considering the case of a static thermal cloud, our results will be correct only for the modes 
where the thermal cloud does not play an important role, i.e., for the out-of-phase modes 



22| , |23|| . In the experiments of Jin et al. |37J, this turns out to be the quadrupole mode. We 
calculate the frequency of the quadrupole mode for this experiment, by means of a fit to the 
experimental data for the damping. 

The frequency and damping of the collective modes are, as measured in experiment, 
averaged quantities. Therefore, we first write down the equations of motion for the averages 
of the gaussian variational parameters and the number of particles in the condensate. The 
equations of motion for the average of the gaussian widths read 

\mNM4M + ^r^ K H = -^ (q(t) ' Nc{t)) - (92) 

For notational convenience we omit the brackets (• • •) denoting the noise average of a stochas- 
tic quantity, and denoted the averages of the gaussian variational parameters simply by qj(t), 
where j equals x, y or z. The average equation in Eq. (|92"D is obtained from the Langevin 
equation for qj(t) by simply leaving out the noise term. This can be done because the noise 
in the Langevin equation does not induce a drift term for the average. This follows directly 
from the Fokker-Planck equation, with the use of partial integration. Since we want to 
describe a perturbation around a static equilibrium, the Keldysh self-energy will be time 
independent to a good approximation, and we thus drop its explicit dependence on time. 
For a description of the collective modes, we also have to consider variations in the average 
number of particles of the condensate, caused by the excitation of a mode. This means that 
we also have to consider the rate equation for the average number of particles 

dNJt) /3. nK , M lAr/N iS K 

—£+ = -|zE K [/i c (t) - fi] N c (t) + — . (93) 

In writing down this equation, we again left out the brackets (•••), which denote averaging 
over different realizations of the noise in the stochastic rate equation in Eq. (pq) . The last 
term on the right hand side of Eq. fl9~5|) is a so-called noise-induced, or spurious drift term. It 
arises because in Eq. fl88D we are dealing with multiplicative Stratonovich noise. Without this 
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drift term, the equilibrium number of particles predicted by Eq. ( |93|) would not be correct, 
as we will see lateron. Note that the average of the stochastic rate equation in Eq. ( p3|) is 
very similar to the result obtained by Gardiner et al. |38| . However, their expression for the 
chemical potential of the condensate is different since they do not consider a gaussian ansatz, 
and they also have not made the 'classical' approximation to the fluctuation-dissipation 
theorem. 

To obtain a description of the collective modes of the condensate, we have to linearize the 
equations in Eqs. ( p2|) and fl93|) around their time-independent equilibrium solutions. Let us 
therefore put qj(t) = qy -\- 5qj(t) and N c (t) = N + 8N(t) and substitute this in Eqs. fl92"|) 
and (PBp. Equating the zeroth-order terms after linearization results for the average rate 
equation in 



J dV(qM,N ) ^ 

1 \ dN c », 



1 



iV = 0. 



(94) 



From this equation the equilibrium number of particles can be calculated. It is however much 
more convenient to use the number of particles in the condensate as experimental input, and 
calculate the chemical potential // such that Eq. (|94]) is satisfied. The equilibrium conditions 
for the gaussian variational parameters read 



dV(^°\N ) 



0. 



(95) 



from which q}"> can be calculated. For the noninteracting case, where a = in the potential, 



the above equations result in qj = Jfi/mojj as expected. In this case N is given by 



iV n 



l3(-(u} x + u y + u z ) -ft) 



-1 



(96) 



which is the correct equilibrium ground state occupation number of a non-interacting Bose 
gas, within the 'classical' approximation. Note that without the noise-induced drift term in 
the rate equation for the average number of particles the correct equilibrium would not have 
been obtained. 

The linearized equations of motion for the deviations are found by equating the first 
order contribution on the left and right-hand side of Eqs. ( |92"D and (J93|) after linearization. 
The equation for the deviation in the equilibrium number N of condensate atoms due to 
the excitation of a collective mode is given by 



5N(t) = -r 5N(t) -J2 a j S Qj(t)- 



(97) 



Here, the parameter T is given by 



P 



r = -^e k 

2 



1 d 2 V(q(°\N ) 

+ JV - 



PN 



dm 



(98) 



where we eliminated the chemical potential \x by using Eq. fl94|). Physically, T describes the 
lack of detailed balance between the thermal cloud and the condensate due to an excitation 
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of a collective mode of the condensate. In general, this lack of detailed balance will cause 
damping, and will alter the frequency with respect to the undamped case. The parameters 
atj are given by 

and represent the response of the fluctuations in the number of particles due to a deformation 
of the condensate in the j-th direction, and vice-versa. 

The linearized equations of motion for the deviation of the gaussian variational param- 
eters take the form of damped harmonic equations 

8q 3 (t) + VjSqjit) = -a 3 /a 6N(t) - fij%(t) - ]T n%Sq k (t). (100) 



The damping rates T 3 are given by 

^ = - rmz? W> ( 101 ) 



2m(qf^ 2 



and the frequencies fL and flj k read 



2 2 d 2 V( q ^,N ) 2 2 d 2 V(^,N ) 

j mN 8q] ' jk mN dq 3 dq k ' l ' 

In Eq. ( [LOOp we introduced the parameter a = mNQJ3iTj K j '4 for later convenience. Physically, 
the damping rates Tj arise because of collisions between thermal atoms and condensate 
atoms. This causes damping of the collective modes on the condensate, and also alters the 
frequencies with respect to the results obtained without damping. It should be noted here, 
that all the parameters mentioned above can be calculated microscopically, by using the 
expression for the Keldysh self-energy given in Eq. (|62}). 

To obtain the eigenmodes of Eqs. (|97|) and ( [LOOP we have to consider solutions of the 
form (SN(t), <5q(£)) = (SN(0), 5q(0))e~ tujt . For such solutions we can rewrite these equations 



as 
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0. (103) 



This matrix equation only has non-trivial solutions if the determinant of the above matrix 
is equal to zero. Solving this condition for uj results in general in complex frequencies 
oo — o>r — iuj\. The real part cur then gives the frequency of a collective mode, whereas the 
imaginary part uj\ gives the damping of this mode. 
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1. Isotropic trapping potential 

We consider now the case of an isotropic trapping potential \^ ext (x) = |mco>QX 2 for 
a discussion of the frequencies and damping rate of the low-lying collective modes of the 
condensate. Because of the spherical symmetry of the condensate in equilibrium, the number 
of parameters reduces and the eigenvalue equation simplifies significantly. For the frequencies 
and damping rates we have: 
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T x 


= r, 


= T Z 
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(104) 

For the parameters atj we have the same simplification, and we denote these parameters 
by a r . With these simplifications, the eigenmodes can be calculated analytically. One 
mode is doubly degenerate, with eigenvectors (0, 1, —1, 0) and (0, 1, 0, —1), and is thus the 
quadrupole mode. From the form of the eigenvectors of the quadrupole mode, we are led to 
the important conclusion that the number of particles in the condensate is constant for this 
mode. Physically this can be understood from the fact that the motion of the condensate 
is 'volume-preserving' in this case: as one direction shrinks the other one expands and vice 
versa. This motion does not lead to a change in the chemical potential of the condensate, at 
least to a linear approximation, and therefore does not affect the average number of atoms 
in the condensate. The complex frequency of the quadrupole mode is given by 



^uad = yn? - n 2 „ - (y) -iir P . (105) 

Note that from this expression it is clearly seen that the damping also affects the frequencies 
of the collective modes [^3[ . 



The frequency and eigenvector of the monopole mode can also be calculated analytically 
for the isotropic case. However, because of the rather formidable expressions involved, we 
omit them here. For the monopole mode the number of atoms in the condensate is not 
constant, but oscillates out of phase with the spatial degrees of freedom of the condensate. 
This can be understood from the fact that the monopole motion leads to a global increase 
in the density of the condensate, and therefore affects the detailed balance of the condensate 
with the thermal cloud. In the case where we ignore the fluctuations of the number of atoms 
in the condensate, and take a = a r = T = 0, the expression for the complex frequency of 
the monopole mode is given by 



Wmono = \jW r + 2QI - (y) - Uv r . (106) 

Comparing the results in Eqs. ( |105| ) and ( |106| ) we see that the damping rate of both modes 
is equal in first order in T r , at least within our variational approximation. 

We now turn to an explicit calculation of the frequencies and damping rates of the 
quadrupole and monopole mode in an isotropic trap. We have used the same parameters as 
Williams and Griffin [ PB] , P6] , and thus have taken the trapping potential frequency equal to 
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uo/2tt = 10 Hz. The calculations are performed for 87 Rb, which has a scattering length of 
a = 5.7 nm. We take the total number of atoms equal to N tota i = 2 x 10 6 . Since the number 
of atoms in the condensate is large at most temperatures below the critical temperature, we 
are mainly in the Thomas- Fermi regime, and can neglect the kinetic energy of the condensate 
atoms with respect to their mean-field interaction. Therefore, we have used a Thomas- Fermi 
profile for the condensate to calculate the collision integral in the expression for the Keldysh 
self-energy in Eq. (|B2|). As a result the Keldysh self-energy turns out not to be constant over 
the size of the condensate in this limit, contrary to our assumption in the derivation of the 
stochastic equations for the variational parameters. To compensate for this effect, we have 
calculated /lE K by taking a volume average of fr£ K (x) over the size of the condensate. We 
report our results for the number of condensate atoms, the frequencies and the damping rates 
as a function of the reduced temperature T/Tbec, where Tbec is the critical temperature 
for an ideal Bose gas. 

The procedure for calculating the number of atoms in the condensate as a function of 
temperature is as follows. For a given number of condensate atoms N we first calculate 
the average condensate density profile and the chemical potential from the time-independent 
Gross- P it aevskii equation. In the Thomas- Fermi limit the average condensate density profile 
is given by 

K0(x))| 2 = ^(/i-^ xt (x)), (107) 

with a condensate chemical potential 

^ (15iV a/Z) 2/5 , (108) 



where I = Jh/mtoo is the harmonic oscillator length. Clearly, the condensate density can 
not be negative, so Eq. ( |107|) is only valid if the condensate density is positive, otherwise 
it should be taken equal to zero. At |x| = iJ fi/tuvo = -Rtf the condensate density is equal 
to zero. Next we calculate the number of atoms in the thermal cloud with the value of the 
chemical potential determined by Eq. fll08| ) using 



A thermal = J dx J -^LjV( e (x, k)), (109) 

with A(e(x, k)) given by Eqs. fl5"T|) and (pB|). We repeat these steps for a variable number 
of condensate atoms until N + A thermal = iVtotai- The result of the calculation is shown in 
Fig. m together with the result for an ideal Bose gas in the thermodynamic limit. Using 
this result for N as a function of the temperature T, we have calculated the Keldysh self- 
energy as a function of temperature using the expression in Eq. ( |62|) ||39|1 . In Fig. H the 
function |/3fr£ K (x)| is shown, for T = 0.5Tbec- It is clear from this figure that the Keldysh 
self-energy is not constant over the size of the condensate, and even diverges in the Thomas- 
Fermi approximation at |x| = -Rtf- The equilibrium values of the variational parameters 



were calculated using Eq. (95). Subsequently, the various parameters were calculated from 



Eqs. fl9"S|), (|9"9"D, ( |101| ) and ( |102p . The complex frequency of the quadrupole was calculated 



from Eq. ( |105| ), and the complex frequency of the monopole mode was calculated from the 
corresponding analytical expression, which we have omitted here. 
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The results for the frequencies are presented in Fig. |3|. The dashed lines are the T = 
frequencies obtained by Stringari flOf , in the Thomas- Fermi limit. Since the calculations 



are in the collisionless limit, where |/%X K | <C 1, the results are essentially T = results 
for a variable number of condensate atoms. In the Thomas-Fermi limit these frequencies 
are independent of the number of atoms in the condensate, but below iVo ~ 10 4 condensate 
atoms, where the Thomas- Fermi approximation starts to break down, the frequencies deviate 
from the T = results, as is seen in Fig. ^. Although for the monopole mode the full 
expression for the frequency and damping involves the parameters T, a r and a, which are 
related to the fluctuations in the number of condensate atoms due to the excitation of a 
collective mode, this hardly affects the results for the frequencies and damping of this mode. 
We therefore conclude that the fluctuations in the number of condensate atoms during the 
excitation of a collective mode hardly affect the damping and frequency of this mode, and 
that the expression in Eq. (|106|) is valid as long as |/3frS K | <C 1. 

The results for the damping rate are shown in Fig. [|, together with the result for |/3^E K |. 
Within our variational approximation, the damping rates for both the quadrupole and the 
monopole modes are found to be the same in the collisionless regime considered here. As 
clearly seen from Fig. |], the damping rate increases with increasing temperature. This is 
because the density of the thermal cloud becomes larger with increasing temperature, and 
there are therefore more collisions between condensate and thermal atoms, which cause the 
damping. Williams and Griffin |35 j have also calculated the damping of the monopole mode 



for a condensate in a spherical trap in presence of a static thermal cloud. These authors 
have generalized the wave equation derived by Stringari ^(| to nonzero temperatures, simi- 



lar to our result in Eq. (|T2|), albeit that their work does not obey the fluctuation-dissipation 
theorem as mentioned previously. They have calculated the damping of the monopole mode 
in perturbation theory, considering the damping as a perturbation parameter. Our results 
for the damping of the monopole mode have the same order of magnitude as their results. 
There are however some qualitative differences. We find that the damping rate increases 
very slowly with temperature for a large temperature regime T < 0.95Tbec 5 and then in- 
creases dramatically as the temperature approaches the critical temperature. Williams and 
Griffin find that the damping rate increases much more gradually with increasing temper- 
ature. These differences are probably due to the fact that these authors take into account 
that the collision integral in Eq. fl62|) has a position dependence. In Ref. |36| Williams and 



Griffin improve upon their Thomas-Fermi calculation for the damping rate by using the 
Bogoliubov-deGennes equations that follow from linearizing the Gross-Pit aevskii equation 
with an imaginary term. In this case, their results for the damping also show a dramatic 
increase as the temperature reaches the critical temperature and the Thomas-Fermi ap- 
proximation breaks down. In this latter work it is also found that the damping for the 
quadrupole mode and the monopole mode are slightly different. The fact that we find that 
the damping for both the monopole mode and the quadrupole mode are equal to first order 
in T r is a result of neglect of the spatial dependence of the Keldysh self-energy Williams 
and Griffin also find that there is no first-order correction in the damping to the real part 
of the frequencies, a conclusion consistent with Eqs. ( |105|) and ( |106|) . 

Summarizing, we have calculated the damping and frequencies for the quadrupole mode 
and the monopole mode of a condensate in a spherical trap. Our results differ slightly both 
qualitatively and quantitatively from the theoretical results found by Williams and Griffin 
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|, [3q| . These differences are probably mostly due to the fact that the calculations are 
performed for a large number of atoms in the trap, which implies that the collision integral 
involved in the calculation has a significant position dependence, which our variational 
approach does not properly account for. However, for a smaller amount of atoms in the 
trap, we believe that our method should give accurate results, and goes in principle beyond 



the perturbation theory considered in Refs. [p5| , [36 



2. Anisotropic trapping potential 

We now calculate the frequency of the m = 2 quadrupole mode, where m is the azimuthal 
quantum number of the angular momentum, for the experimental parameters of Jin et al. 
| |37| |. In this experiment, one loads 87 Rb atoms into an anisotropic trap, with radial frequency 
uj r /2iT = 129 Hz, and axial frequency uj z /2-k = 365 Hz. Although the equilibrium shape of 
the condensate is now anisotropic, the expression for the frequency of the quadrupole mode 
found in the isotropic case in Eq. ( |105|) , turns out to be also correct for the m = 2 quadrupole 



mode. The parameters fl r , fl rr , and T r , are now given by 



i L rr r= i L X y ^^yxi 



T r = Vt x = Vt yi (110) 

which follow from the axial symmetry of the condensate in equilibrium. It follows from the 
expression for the complex frequency of the quadrupole mode in Eq. (|1(J5|) that the complex 



frequencies lie on a circle of radius |co>g Ua( i| = JVl% — f^ r ~ \/2uj r . To test the validity of our 
expression for the frequency of the quadrupole mode, we have plotted the experimental data 
points taken from Ref. |37j in the complex cu-plane. In Fig. [| the result is shown, together 
with a circle of radius \/2uj t . The good quantitative agreement with experiment, clearly 
visible in Fig. ||, implies that the expression in Eq. ( |105| ) for the frequency and damping 
of the quadrupole mode is correct, even in the hydrodynamic regime, where |/3fr£ K | ^> 1. 
This may in first instance come as a surprise, since our variational approximation to the 
stochastic non-linear Schrodinger Eq. (|60"D was derived in the collisionless regime, where 
|/%£ K | <C 1. Apparently, the relation (wquadl — V2u r is also valid in the hydrodynamic 
regime. This can be understood from the fact that this relation for the complex frequency is 
quite general for a damped harmonic oscillator and we expect on general grounds that the 
quadrupole mode of the condensate can be described in this way, both in the collisionless 
and in the hydrodynamic regime. 

To determine |/3frS K | as a function of temperature, we have fitted the imaginary part of 
Eq. ( |105| ) to the experimental data for the damping presented in Ref. ||37|| . From this fit, we 
have calculated the dimensionless parameter |/3fr£ K |, using Eqs. (|101|) , (|102| ), and (|105|) . The 



results of this fit are presented in Fig. || The value for |/3fr£ K | found in this manner is then 
used to calculate the real part of Eq. ( |105| ), i.e., the frequency of the collective mode. In both 
the fitting of the damping, and the calculation of the frequency, we used the fit presented 



in Ref. [23| to determine the experimental values for the number of condensate atoms as a 
function of the temperature. The result for the frequency is presented in Fig. [7], together with 
the frequency calculated from Eq. ( |105|) , with T r = 0, i.e., the zero-temperature frequency 
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for a variable number of condensate atoms. The experimental points are also shown. In 
Fig. [7| a good quantitative agreement with the experimental results is found. Fig. |7] also 
shows clearly that at nonzero temperatures the damping seriously affects the frequency of 
the quadrupole mode, as also found by Al Khawaja and Stoof [p3fl . However, a microscopic 
calculation of fr£ K in the Thomas- Fermi limit for the experimental conditions of interest, 
by means of a volume average of Eq. (J3^) over the size of the condensate, turns out to 
give values for |/5/iS K | which are approximately one order of magnitude too small to explain 
the experimental data for the quadrupole mode. There are several possible reasons for this 
discrepancy. One possible reason is that the calculation of KL K as an average over the 
size of the condensate of fr£ K (x) is not a good approximation. Since fr£(x) becomes large 
compared to fr£ K (0) at the edges of the condensate, where the depletion of the thermal cloud 
due to the condensate's mean field is relatively small, the position dependence of fr£ K (x) 
is of importance, in particular since the density fluctuations occur for the quadrupole mode 
precisely near the edges of the condensate. This is clearly shown in Fig. |2|, which shows that 
the Keldysh self-energy even diverges at |x| = R^p in the Thomas- Fermi limit. With respect 
to this remark, we refer to future work concerning the full numerical solution of the average of 
the Langevin equation in Eq. (0), to investigate the importance of the position dependence 
of fr£ K (x), and a comparison of this numerically exact approach to the variational method 
developed here. Another possible reason for the discrepancy with the experimental results, 
is the presence of other sources of damping, such as Landau and Beliaev damping, which 
have not been included in our calculations. 



B. Condensate growth and collapse 

Although so far we have focused on repulsive interactions, and thus a positive scattering 
length a, we consider in this section the case where the scattering length is negative. In this 
case the condensate energy in the gaussian approximation becomes 

v^)^{— ]+xm N^--n-^-. (in) 

From this potential it can easily be seen that there is only a metastable condensate possible 
if the number of atoms in the condensate is smaller that a certain critical value. This 
is illustrated in Fig. || where we show the potential in Eq. (ITT), for several values of the 



number of atoms in an isotropic condensate. For an isotropic trap, the maximum condensate 
number iV max turns out to be given by the condition 

iVmax^ < ^ * 0.67. (112) 

If the number of atoms in the condensate is above this value, the potential in Eq. (|111|) 
has no (meta) stable minima. If the number of atoms is smaller than iV max , the potential 
has a metastable minimum, and the condensate can start to collapse by overcoming the 
metastable energy barrier by either macroscopic quantum tunneling or thermal fluctuations 
T3| . The stability condition for the condensate, found by a full numerical solution of the 
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Gross-Pitaevskii, turns out to be N c a/l < 0.58 [41], and thus we see that the gaussian 
approximation is only 16 % off. 

The first experiments on Bose-Einstein condensation in a gas with attractive interactions 
where performed with 7 Li [[14]], which has a negative scattering length of a ~ —1.45 nm. 



In these experiments, the gas is evaporatively cooled below the critical temperature, which 
causes the condensate to undergo several growth and collapse cycles before relaxing to a 
metastable equilibrium |T5 ]. Because of the stochastic initiation of the collapse, one could 



not make a sequence of destructive measurements. However, a statistical analysis revealed 
that during a collapse the number of condensate atoms is not reduced to zero, but that the 
collapse stops, presumably because of elastic and inelastic collisional loss processes, when 
the number of atoms in the condensate is about 200 ||16|| . In a recent experiment, one was 
able to make a sequence of destructive measurements by (partially) dumping the condensate 



by a two-photon pulse ||17|| , and thus start each measurement with approximately the same 
initial number of condensate atoms. 

Using the Langevin equations for qj(t), and the stochastic rate equation for the number 
of atoms in the condensate, we are able to describe this experiment. To do so in the 
easiest way, we want to model the thermal cloud by an equilibrium Bose distribution of a 
noninteracting gas. However, numerical solutions of the quantum Boltzmann equation for 
these experiments have shown that the thermal cloud is not equilibrium, but can be well 
modeled by a distribution function given by |16 



exp [/?(// - fj)] _ A 

/U exp[/3(e -fi)]- 1 exp[/3(e-/i)]-l' l ' 

At high energies, / has the form of a Boltzmann distribution with chemical potential \J and 
temperature l/ksP- At low energies / has the form /(e) = A//3(e — fi), which is precisely 
the low-energy tail of a Bose-distribution with effective temperature A/k^fi. Therefore, we 
conclude that Eq. fl54"|) is too a good approximation still valid for the distribution function 
given by Eq. (|113|) , if we make the replacement (3 — > (3/ A in the operator on the right hand 
side of Eq. (|5"4]). In this manner we obey the fluctuation-dissipation theorem, and have also 
accounted for the fact that the distribution function of the thermal cloud is a nonequilibrium 
distribution function. We expect that this approximation will give quantitatively correct 
results for the condensate growth rate, since it is particularly good for the low-lying energy 
levels, which dominate the condensate growth. 

Before comparing to the experimental data reported in Ref . |DJ , we discuss some aspects 



of the numerical solutions of the stochastic rate equation coupled to the Langevin equation, 
for a condensate which has initially no atoms. In the experiment performed by Gerton et 
al. this would correspond to the situation where the condensate is dumped completely. The 
stochastic rate equation in Eq. flggp is well suited for this purpose, since it also contains 
fluctuations, which initiate the growth in this case. Without these fluctuations, the growth 
rate of the condensate would never become nonzero. These initial conditions for condensate 



growth are different from the experiments conducted by Miesner et al. [l"Hfl, in which the 
condensate growth is observed after evaporatively cooling the gas. In this case, the ground 
state already has a rather large nonzero occupation number above the critical temperature, 
which causes a growth process dominated by Bose-stimulation. Therefore, for a theoretical 
description of the condensate growth it is not so essential to include fluctuations in this case 
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43| , |44|| . We perform our simulations for the experimental conditions reported in Ref. [17 . 



The trap frequencies are given by uj r /2iT = 151 Hz and uj z /2ir = 131.5 Hz. We consider a 
thermal cloud with 70000 atoms at a temperature T = 170 nK. The parameter A is taken 
equal to 4. These values correspond to typical experimental conditions |4l|. The Keldysh 



self-energy fr£ K is calculated using Eq. (|62|), with the nonequilibrium distribution function 
/(e) given by Eq. ( |113| ). The mean- field effects of the condensate on the thermal cloud 



are neglected, an approximation which will certainly be valid in the initial stage of the 
condensate growth, when the condensate is small. Moreover, we take the chemical potential 
of the condensate (/x c (x, £)) in Eq. ( |5^D equal to zero. Since the density of thermal atom will 
be the largest in the center of the trap, and the condensate is small in this case, we do not 
perform an average to calculate fr£ K , but simply take hH K = fr£ K (0). 

We solve the Langevin equations coupled to the stochastic rate equation, using standard 



numerical techniques for stochastic differential equations ||46|| . Since the initial number of 
condensate atoms is equal to zero, we put at time t = the values of the gaussian variational 
parameters equal to qj = Jh/mujj, which is their equilibrium value in the limit where the 
number of condensate atoms approaches zero. A slight subtlety arises in the use of the 
stochastic rate equation. One has to realize that the chemical potential of the thermal 
cloud in this equation is measured with respect to the energy of the lowest excited level, 
since we want the distribution function to describe the non-condensed atoms only. This 
means that we should use in the rate equation the chemical potential found in matching 
the distribution function in Eq. (|113|) to the number of thermal atoms, and add the energy 
of the lowest excited level to it. This is immediately clear when we write the number of 
atoms in the thermal cloud as a sum over occupation numbers, instead of an integral over 
energy. Fig. |9] shows the results of our simulations. In Fig. |9] (a), we plot the number of 
condensate atoms as a function of time, for the solutions of the Langevin equations and 
the stochastic rate equation for three different realizations of the noise. We assume that 
during the growth and subsequent collapse of the condensate, the distribution function of 
the thermal cloud is not affected. The maximum number of atoms in the condensate is for 
the parameters under consideration here equal to iV max ~ 1470 atoms, within the gaussian 
approximation. This means that during one growth-collapse cycle the number of atoms in 
the thermal cloud is reduced by only approximately 2 %, and therefore this approximation 
seems valid for the description of one growth-collapse curve. Once a collapse is initiated 
by the noise in the Eqs. flSU] ) and (|55|), we model the collapse by putting the number of 
condensate atoms instantaneously equal to a gaussian random number with a mean value of 
200 and a deviation of 40, which corresponds to the 20 % systematic uncertainty reported 
in |TJ|, and the variational parameters qj(t) equal to their corresponding equilibrium values, 
given by Eq. (|95"|). The collapse occurs on a time-scale 0(l/u), which is much faster than the 
time-scale on which the condensate grows. Since we are interested in the growth process here, 
and not in the loss process which stops the collapse, this appears a reasonable way to model 
the collapse. Fig. [| (a) clearly shows that when the number of condensate atoms approaches 
the maximum number iV max , the condensate tends to collapse. We found that the collapse is 
initiated stochastically by the fluctuations in the number of atoms in the condensate, which 
cause density fluctuations that cause the condensate to overcome the macroscopic energy 
barrier and start the collapse. Since our description only includes thermal fluctuations, 
we may ask if macroscopic quantum tunneling might be of importance. However, previous 
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work has shown that decay by thermal fluctuations is the main decay mechanism for the 
temperatures of interest | 15| . In Fig. || (b) we plot the number of atoms in the condensate 



as a function of time, averaged over different realizations of the noise. The red line is an 
average over 5 different realizations of the noise. A growth-collapse signature is still visible in 
this curve, although the stochastic growth process and initiation of the collapse has led to a 
dephasing of the moment of collapse. The green line shows an average over 10 realizations of 
the noise. Although the initial growth is clearly visible in this curve, the collapse can hardly 
be seen from this average, since the noise has led to an almost complete dephasing, and the 
collapse is 'averaged out'. Finally, the blue curve shows an average over 1000 realizations of 
the noise. No signature of the collapse is visible in this curve, because the averaging leads 
to a complete dephasing of the moment of the collapse. 

We now discuss the simulation of the experiments performed by Gerton et al. JT7J. To 



make a comparison with experiment, one has to realize that each data point is obtained as 
an average over 5 or 10 individual runs. Since the condensate number is probed by means 
of a destructive measurement, each experimental curve should not be viewed as an average 
of curves. Instead, each point is an average over different experimental runs, and the time- 
correlation between different experimental points is only caused by the initial conditions, 
which are approximately the same for each experimental run. To simulate this experiment by 
means of a numerical solution of the Langevin equations in Eqs. fl86|) and (p8|), we therefore 
have to let the numerical solutions evolve up to a certain point in time, and then make a 
numerical measurement. We then average over 5 or 10 measurements to obtain a data-point 
and its uncertainty, and repeat this procedure at a different measurement time. In this 
way, we are certain that each individual solution of our stochastic equations, is done for a 
different realization of the noise. Note that this procedure is very reminiscent of the method 
of Monte Carlo simulation. 

We have done simulations for the three different experimental situations presented in 
Ref. p!7| . The results of our simulations are presented in Fig. 110, as red triangles. The 
experimental data-points are also shown, and denoted by black circles. The Keldysh self- 
energy was calculated as in the simulations described above, using the averages of the full 



experimental data on the number of atoms in the thermal cloud and their temperature |H5 



For the parameter A of the nonequilibrium distribution function in Eq. ( |113| ) we use the 



average of the fits obtained by Gerton et al. [fJ5[] . For Fig. |10] (a) and (c) this corresponds 
to a thermal cloud of approximately 65000 atoms at a temperature of T = 170 nK. The 
parameter of the nonequilibrium distribution is equal to A = 4 for these simulations. For 



Fig. [10| (b) the thermal cloud contains approximately 100000 atoms at a temperature of 
T = 200 nK. The parameter A = 2 in this case. For Fig. |TU] (a) and (b) the averages are 
taken over 5 runs for each data-point, whereas for (c) 10 runs are used. The error bars 
denote the uncertainty in the average. In Fig. [H] (a) and (b) the initial number of atoms 
was take equal to iV c (0) = 100, and iV c (0) = 438 respectively. For Fig. [10] (c) we have 
taken iV c (0) = 0, since the condensate was in this case dumped completely to within the 
experimental resolution ||17||. 



The results of our simulations presented in Fig. |10| shows good agreement with experiment 
for the initial stage of the growth, where the condensate is small. In particular, Fig. |1(] (a) 
show good agreement in the initial stage where N c < 400 atoms, whereas Fig. pi] (c) shows 
good agreement in the regime where N c < 600 atoms. This is to be expected, since the 
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gaussian ansatz is a very good approximation for a small number of atoms in the condensate, 
whereas it becomes worse for a larger number of atoms in the condensate. The fact that 
the error bars of the experimental data points have the same order of magnitude as the 
error bars on our simulations, indicates that the fluctuations, i.e., the noise in our stochastic 
equation, have indeed the correct order of magnitude. As mentioned in the discussion of 
the individual solutions of our stochastic equations, we find that the collapse is initiated by 
fluctuations in the number of atoms, and that these fluctuations thus lead to a dephasing of 
the moment of the collapse. In principle also the fluctuations in the initial number of atoms 
in the condensate lead to dephasing. However, the uncertainty in iV c (0) is small compared 
to the uncertainty in N c (t) at later times, and we therefore conclude that fluctuations in the 
initial conditions for the condensate are presumably less important for an understanding 
of the dephasing of the moment of the collapse. Moreover, there are also fluctuations in 
the properties of the thermal cloud for each individual experimental run, which we have not 
taken into account. With respect to this point we also note that our method does not display 
the saturation in the growth rate, observed in numerical solutions of the quantum Boltzmann 
equation [17]]. This effect is also observable in the experimental data in Fig. [K] (c), where 
the growth is exponential in the first stage, but then turns linear. This saturation in the 
growth rate is caused by the fact that the condensate mostly grows from the low-lying 
excited states, which in turn have to be feeded by collisions in higher energy states, which 
are not Bose enhanced. Since in our simulations the thermal cloud is taken to be static, 
our simulations do not display this saturation effect. In conclusion, we like to point out 
that, to make a sensible quantitative comparison to the experimental results in the whole 
time-domain, we have to compare the converged averages of both the experimental runs and 
the theoretical simulations. This is because of the fact that the fluctuations are so large, 
that each individual growth curve can differ substantially, as seen from Fig. [5] (a). In turn, 
this leads to averages that can differ qualitatively, depending on the number of runs one 
averages over, as also clearly seen in Fig. |9| (b). 

VI. CONCLUSIONS 

We have presented a Fokker-Planck equation that describes the nonequilibrium dynam- 
ics of an atomic Bose-Einstein condensate. We have discussed an approximation to this 
Fokker-Planck equation, which assumes that the thermal cloud is close to equilibrium. Its 
corresponding Langevin equation has the form of a stochastic nonlinear Schrodinger equation 
with complex gaussian noise. Both the Fokker-Planck equation, and the Langevin equation 
obey the fluctuation-dissipation theorem, which ensures that the condensate relaxes to the 
correct equilibrium. We have also presented the hydrodynamic formulation corresponding 
to this stochastic non-linear Schrodinger equation, in which the condensate is described in 
terms of its density and its phase. These turn out to obey a stochastic continuity equation 
and a stochastic Josephson equation, respectively. To make analytical progress, we have 
then extended the variational calculus, commonly applied to the Gross-Pit aevskii equation, 
also to the case of the stochastic non-linear Schrodinger equation. The equations of motion 
for the variational parameters turn out to be identical to the Langevin equations describing 
the Brownian motion of a particle in a potential. These equations are then coupled to a 
stochastic rate equation for the number of atoms in the condensate. We have applied these 
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equations to calculate the damping and frequencies of the collective modes of the condensate, 
and to obtain a description of the growth-collapse curve of a condensate with attractive in- 
teractions. However, there are much more applications possible with the variational method 
presented here. With a slight extension, it can also be used to calculate the frequency and 
damping of the scissor modes of the condensate |57]J5!| , at nonzero temperatures. Moreover, 
applying the method to a Thomas-Fermi density profile, we can obtain a simple description 
of the growth of a condensate with repulsive interactions. The treatment of the dissipative 
dynamics of vortices and other topological excitations such as skyrmions |^9| , is also feasible 
within this variational method. 

The only quantity that characterizes the thermal cloud is, in our approach, the Keldysh 
self-energy ft,E K (x, t). The parameter that enters the equations of motions for the variational 
parameters, turns out to be some spatial average of this quantity. In our calculations pre- 
sented here, we have taken an average over the size of the condensate, when calculating the 
frequency and damping of the collective modes. Future work will include a numerical solu- 
tion of the stochastic non-linear Schrodinger equation, to investigate the importance of the 
spatial dependence of the Keldysh self-energy, which we have not taken into account here. 
Nevertheless, we believe that the variational method presented here, provides a satisfying 
picture of the nonequilibrium dynamics of a Bose-Einstein condensate at nonzero tempera- 
tures. Moreover, as we have shown, calculations done within this variational approximation, 
lead already to a good agreement with experiments on collective modes and on condensate 
growth. 
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APPENDIX A: AMPLITUDE AND PHASE VARIABLES 

The condensate is often described in terms of density and phase variables, by making 
a canonical transformation <fi = ^/pe* 6 '. In this appendix we discuss the derivation of the 
Langevin equations of motion for p(x, t) and 8{x,t). For simplicity, we first discuss the 
single-mode version of the probability distribution in Eq. ([0]) in the non-interacting case. 
So we consider a probability distribution for a single-mode complex order parameter, which 
reads 



P[0*,0;£] 




d 



ih -Ep + A* - ^c + iR I <f>(f) 



(Al) 



Physically, this probability distribution describes the nonequilibrium dynamics of a non- 
interacting Bose-Einstein condensate, in contact with a thermal cloud characterized by a 
Keldysh self-energy 7i£ K , with inverse temperature (3 and a chemical potential \x. The 

35 



energy per particle in the single-mode system is equal to // c . The dissipation R is again 
related to the Keldysh self-energy by the fluctuation-dissipation theorem in Eq. fl54]) , which 
in this simple case reads 

iR = -^hE K ^ c -ij]. (A2) 

From the previous sections we know that the probability distribution P[(f)*, 0; t] is generated 
by the Langevin equation 

in d^_ = {fic _ fi _ iR) 0(t) + ^ (A3) 

where the complex noise has a time correlation function given by 

(v*(t')v(t)) = ^ K S(t'-t). (A4) 

As explained in the first section, the noise r)(t) has to be interpreted as an Ito process. We 
can again derive the Fokker-Planck equation for P[<f>*,<f>]t\ by noting that it is in fact the 
Schrodinger equation in the position representation. We will do this in some detail once 
more, to make clear the different steps of the derivation. The first step is to determine the 
momenta conjugate to the coordinates and 0*. Since we have a lagrangian equal to 



LWA] 



Un— + /i-/i c + zi?J0(t'; 



we can define the momentum conjugate to in the usual way 



(A5) 



p *=dj = mt[- ih dt-^- zR+fM )' (A6) 

with the complex conjugate expression for p^*. The second step is to derive the hamilto- 
nian. Although it has in principle ordering problems, we overcome these by noting that in 
the path-integral formulation of quantum mechanics we are always dealing with a normal 
ordered hamiltonian. The normal ordered hamiltonian, i.e., with the momentum operators 
positioned left to the coordinate operators, is now given by 

H[p^A;Pr^*} = Pit + Pr<l>* - L [<t>* ,<!>]• ( A7 ) 

The last step towards the Fokker-Planck equation is to quantize the hamiltonian, and to write 
down the Schrodinger equation in the position representation. So we have p^ = —ihd/d<p, 
and similarly p^* = —ihd/d(p*. The Fokker-Planck equation becomes 

^-P[0*, 0; t] = - — (// c - » - iR)<j>P[<j>*, 0; t] 

^-(// c -/i-zi?)0*P[0*,0;£] 



'■" ^S K P[0*,0;t]. (A8) 



2 90*90 
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With the use of the fluctuation-dissipation theorem in Eq. (|A2|) we then show that it has as 
a stationary solution 

P[\<f>\,t -+ oo] oc exp {-PfiHc - nW , (A9) 

which only depends on the amplitude of (f> and <fi*. 

We now repeat the above discussion in terms of amplitude and phase variables, defined 
by yNe l6 . Let us first discuss the equilibrium properties we expect in terms of the number 
of particles N and the phase 9. Since the transformation to N and 9 has a jacobian equal to 
zero, we can just substitute it into Eq. (|A9|), to obtain the equilibrium distribution in terms 
of the number of particles. It is given by 

P[N; t -»• oo] oc exp {-/3(> c - //) A} . (A10) 

We can easily check that the average number of particles in the single mode is in equilibrium 
given 

(7V) - f °° dN P[N;t ^ oo] ~ [fj{fic ^ ■ (AU) 

This is precisely the Bose distribution for /3(/j, c — //) <C 1. If we do not apply the 'classical' 
approximation to the fluctuation-dissipation theorem, as in Eq. (|A2|) , but use the exact 
relation 

lR = _Ift S K (1 + 2 JV( yUc ))- 1 , (A12) 

instead, we find the Bose distribution as the equilibrium number of particles, as expected. 
For the description of a single-mode Bose-Einstein condensate Eq. (|A2|) is in general a good 
approximation, since \i is very close to /x c below the critical temperature. Let us now try to 
derive the Fokker-Planck equation for A and 9. We can do this by substitution of <p — vNe ld 
into the action in the exponent of Eq. (|Al|). Since the jacobian of the transformation is equal 
to one, we simply have 



f N(t)=N ( { -I 

P[N, 9;t]= / d[N]d[9] exp -S[N, 9] , (A13) 

J 9(t)=6 L ft ) 



rN(t)=N ( I 

>8(t) 

with an action equal to 



^•"i = l', dt ' (1^ ("C) + *«■ - tf + ^W) (*« + f mt'))) . (am) 

Naively, we could derive the Fokker-Planck equation from the above path-integral expression 
by going through the same steps as before. If we again apply normal ordering to the 
hamiltonian with respect to N and 9 and their conjugate momenta, the Fokker-Planck 
equation reads 

M \ 2 3/V 2 + h 3N ) l • • ' 

+ {ww + m^-" ) ) p ^^- (A15) 
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This Fokker-Planck equation is however incorrect, since it is easily seen that the equilibrium 
distribution in Eq. ( A10| ) is not a solution of this Fokker-Planck equation. The fact that the 



Fokker-Planck equation in Eq. ( |A15| ) turns out to be incorrect has to do with the fact that 



we have normal ordered the hamiltonian in terms of the variables N and 9. Although normal 
ordering of the hamiltonian Hlp^, 4>]p<i>*, (ft*] did give the correct results, this, however, does 
not imply that we also have to normal order the hamiltonian in terms of N and 9. Let us 
therefore proceed more carefully, and rewrite the Fokker-Planck equation in Eq. ( |A8| ) for 
P[4>*, 4>; t] in terms of N and 9. With the use of the chain rule for differentiation it is easy 
to show that, for a general function / 



Ne 1 "-^ + —= e te -£r, (A16) 

u V dN 2\fN 89 y ' 

with the complex conjugate expression for df/d<p. Substitution of this result in the Fokker- 



Planck equation in Eq. (]A8|) yields the Fokker-Planck equation for P[N,9;t]. It is given 
by 

dt { 2 ON dN + h dN ) [ ' ' J 

+ (ww + m^-ri) p i N >M- (A17) 



Comparison of the Fokker-Planck equations in Eqs. (|A15|) and ( [A 17] ) shows that in Eq. ( |A15|) 



we have misinterpreted the noise on N(t) as an Ito process, whereas in the correct Fokker- 
Planck equation in Eq. ( [A 17] ) we are clearly dealing with a Stratonovich process. Note that 
the same conclusion can also be reached by determining the equation of motion of (N) (t) 
from a variation of the action S[N, 9}. 

From the action in Eq. (|A14|) we can read of the Langevin equations for N and 9. The 
Langevin equation for N(t) is given by 

2R 

T 

(V(t')v(t)) = ^5(t>-t), (A18) 



N(t) = - 2 -^N(t) + 2^m V (t'); 



"(*) . 

N(t) 



4 
and the Langevin equation for 9 reads 

h9(t) = /j, - /j, c 

Ht')v{t)) = l ^8(t'-t). (A19) 

From the above discussion we thus conclude that we have to interpret the noise in the 
Langevin equation in Eq. ( |A18| ) for the number of particles N(t) as a Stratonovich process, 
to achieve the correct equilibrium distribution. This is the main conclusion of this appendix. 
It is straightforward to show that the above discussion generalizes to the case of a multi- 
mode description of the condensate in terms of the complex field 0(x, t). When we make 
the transformation to density and phase variables by setting (ft = y/pe ie , we again have to 
be careful, and interpret the multiplicative noise that enters the equation of motion for the 
density p(x, t) as a Stratonovich process. 
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FIGURES 

FIG. 1. The solid line gives the condensate fraction for a 87 Rb gas of 2 x 10 6 atoms in an 
isotropic trap with ujq/2tt = 10 Hz. The dashed line corresponds to the ideal Bose gas in the 
thermodynamic limit. 

FIG. 2. The dimensionless quantity |/3?iX K (x)|. At |x| = i?TF there is a divergence, when the 
collision integral in Eq. (^) is calculated in the Thomas- Fermi limit. This is indicated by the 
dashed line. The calculation is performed for T = 0.5Tbec with the same parameters as in Fig. [I] 

FIG. 3. Frequencies for both the monopole and quadrupole mode as a function of the tem- 



perature. The dashed lines indicate the zero-temperature results found by Stringari [40 1 , i.e., 
^quad = V%uo an d w mon o = V^uiq- The plot starts at T = 0.8Tbec since for smaller temperatures 
the deviation from these values is neglegible. The parameters are the same as in Fig. |J. 

FIG. 4. The damping rate for both the quadrupole and monopole mode for a condensate in an 
isotropic trap, which is the same for small damping in our variational approximation. The inset 
shows the dimensionless parameter |/37iE |. The parameters are the same as in Fig. |l|. 



FIG. 5. The complex w-plane. In our expression for the frequency and damping rate of the 
m = 2 quadrupole mode in Eq. (105) the complex frequencies lie on a circle of radius \J2uj r . The 
experimental points taken from Ref. [37] are also shown. 



FIG. 6. Fit to the damping rate of the quadrupole mode, as measured by Jin et al. [37|. The 



solid line shows the fit and the experimental points are taken from Ref. [37|. The inset shows the 
value of \/3%T, K \, as calculated from this fit with Eqs. (flOlb , (fof ) and ( |J05l) . 



FIG. 7. The frequency of the quadrupole mode as a function of temperature, calculated with 
Eq. (105), by using the fit shown in Fig. ^. The dashed line show the T = frequency for a variable 



number of condensate atoms. The experimental points taken from Jin et al. |37j are also shown 



FIG. 8. The condensate energy as a function of the condensate width in the gaussian approx- 
imation for three different values of the condensate atom number. If the number of atoms in the 
condensate is larger than -/V max the condensate is unstable, otherwise a metastable condensate is 
possible. 

FIG. 9. (a) Growth-collapse curves of a 7 Li condensate, and (b) their averages. The colored 
lines in (a) display the number of condensate atoms for solutions of the Langevin equation for qj(t), 
coupled to the stochastic rate equation for N c (t), for three different realizations of the noise. In 
(b) the red line correspond to an average over 5 realizations, the green line to 10, and the blue line 
to an average over 1000 different realizations of the noise. The simulations are done for a thermal 
cloud of 70000 atoms with T = 170 nK, and A = 4. 
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FIG. 10. Simulations of the experiment performed by Gerton et al. [17]. The results of the 
simulations are denoted by red triangles, the experimental data is shown as black circles. In (a) 
and (b) each data-point of the simulations is an average over 5 runs, as in the experiment. For (c) 
10 runs per point were done. The error bars in both the experimental data, and the data obtained 
by the simulations, denote the uncertainty in the mean. 
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